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Abstract 

Despite the relative recency of its inception, the theory of compres- 
sive samphng (aka compressed sensing) (CS) has already revolutionized 
multiple areas of applied sciences, a particularly important instance of 
which is medical imaging. Specifically, the theory appears to provide 
an answer to the important problem of optimal sampling in MRI, with 
an ever-increasing body of works reporting stable and accurate recon- 
struction of MRI scans from the number of spectral measurements which 
would have been deemed unacceptably small as recently as five years ago. 
Reducing the number of MR measurements per scan comes to address 
one of the most critical impediments intrinsic in MRI, which is the rel- 
atively slow speed of image acquisition. Although very significant, such 
an improvement may still be insufficient in the cases when a repetitive 
acquisition of MRI scans pertaining to the same volume of interest is re- 
quired. Acquisitions of this type are prevalent in diffusion MRI, in which 
an independent MRI scan is required to encode the strength of water 
diffusion along a predefined spatial direction. Thus, for example, an accu- 
rate delineation of multimodal diffusion profiles by means of high angular 
resolution diffusion imaging (HARDI) requires using between 60 and 100 
(iifl:usion-encoding gradients. This, in turn, is translated into relatively 
long acquisition times, which adversely affects the applicability of HARDI 
for clinical diagnosing. To overcome this limitation, the present paper in- 
troduces a method for substantial reduction of the number of diffusion 
encoding gradients required for reliable reconstruction of HARDI signals. 
The method exploits the theory of CS, which establishes conditions on 
which a signal of interest can be recovered from its under-sampled mea- 
surements, provided that the signal admits a sparse representation in the 
domain of a linear transform. In the case at hand, the latter is defined to 
be spherical ridgelet transformation, which excels in sparsifying HARDI 
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signals. What makes the resulting reconstruction procedure even more ac- 
curate is a combination of the sparsity constraints in the diffusion domain 
with additional constraints imposed on the estimated diffusion field in the 
spatial domain. Accordingly, the present paper describes a novel way to 
combine the diffusion- and spatial-domain constraints to achieve a maxi- 
mal reduction in the number of diffusion measurements, while sacrificing 
little in terms of reconstruction accuracy. Finally, details are provided on 
a particularly efficient numerical scheme which can be used to solve the 
aforementioned reconstruction problem by means of standard and read- 
ily available estimation tools. The paper is concluded with experimental 
results which support the practical value of the proposed reconstruction 
methodology. 

1 Introduction 

The human brain consists of about 10^^ nerve cells that can be subdivided 
into about 1000 different cell types, a complexity that far exceeds that of other 
organs of the body. A further complexity is evident in the way in which the 
component cells of the brain interconnect and function fT . In contrast to other 
types of the cells, each neuron communicates with many target cells by means 
of its distinctive protoplasmatic protrusion, called an axon. Axons with similar 
destinations, in turn, tend to form bundles - known as neural fibre tracts - 
which play a pivotal role in the determination of brain connectivity. Through 
reconstructing the pattern of connectivity of the neural tracts in both healthy 
and diseased subjects, it is therefore possible to obtain an abundance of valuable 
diagnostic information that could be used for early diagnostics of brain-related 
disorders, for assessing the damage caused to the brain by stroke, tumours or 
injuries, as well as for planning and monitoring of neurosurgeries [2j. 

Central to MRI is the notion of contrast, which is typically defined by the 
biochemical composition of interrogated tissue as well as by the morphology of 
its associated parenchyma. Prevalent in MRI practice are the contrasts deter- 
mined by the T1/T2 relaxation times and proton density (PD). Despite their 
exceptional importance to clinical diagnosis, none of the above contrast mecha- 
nisms has demonstrated effectiveness in delineating the morphological structure 
of the white matter. It is only with the advent of diffusion MRI (dMRI) that 
scientists have been able to perform quantitative measurements of the diffu- 
sivity of white matter, based on which its structural delineation has become 
possible [2j-^. Although over the last two decades dMRI has developed into 
an established technique with a great impact on health care and neurosciences, 
like any other MRI technique it remains subject to artifacts and pitfalls [s]. 
While many of such artifacts can be overcome by using advanced hardware de- 
signs and/or more sophisticated imaging protocols [2] Ch.2], one particularly 
critical limitation of dMRI stems from the physics of the acquisition of diffu- 
sion MR images, and therefore is impossible to resolve by operational means. 
Specifically, since collecting the dMRI data requires a repetitive acquisition of 
MR responses from the same volume for a number of diffusion-encoding gradi- 
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ents, it is the relatively long acquisition times that greatly impair the practical 
value of this important imaging modality. Particularly, longer acquisition times 
entail a higher probability for the patient to exercise involuntary motion (typi- 
cally caused by fatigue and/or stress related tremors, swallowing, uncontrollable 
sighing or coughing), which severely affects the quality of dMRI data. Moreover, 
since during the whole duration of the scan the patients are required to remain 
motionless, it is currently deemed ineffective to apply dMRI-based diagnosis 
in paediatrics as well as to patients with dementia or post-traumatic injuries, 
where non-compliance is typical. The problem of long acquisition times also 
hampers the application of dMRI for intra-operational imaging, where it could 
be an irreplaceable tool to use for neurosurgical planning and decision-making 



support 10 12 . Lastly, relatively long scanning times required by dMRI aggra- 
vate the problem of accessibility to MRI equipment. All the above arguments 
suggest that the practical value of dMRI could be improved by shortening the 
scanning times required for acquisition of dMRI data. A particular method to 
achieve such an improvement is detailed in this paper. 

In this work, we adopt a general diffusion model in which each voxel within a 
region of interest (ROI) is allowed to support more than one fibre tract. In this 



case, under some general assumptions (see, e.g., 13 Sect. 3.1] for more details), 
the diffusion signal s(u | r) originating from a voxel with spatial coordinate r G 
containing M(r) fibres can be modelled as 
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M(r) 

s(u|r) = so(r) Q!i(r) exp {-6 (u^A(r) u)} , (1) 



i=l 

where u denotes the spherical coordinate, i.e. 

ue§2 |^^]g3 1 ii^ll^ ^ (2) 

and ai{r) > are positive weights obeying J2i=i'^ cti(i") 1- In ([ll, so denotes 
the diffusion signal obtained in the absence of diffusion encoding (i.e. the so- 
called "bO image"), b is defined as a function of the shape and amplitude of 
diffusion-encoding gradients 15 Eq. 3.18], and {^'^(r)}^^^'"^ are 3x3 diffusion 



tensors associated with the M(r) neural fibres passing through the coordinate r. 
In practical settings, the spherical coordinate u is sampled at K distinct points 
{^k}k=i O'^s'' the unit sphere. In this case, for each u^, MR measurements are 
acquired in the form of a diffusion-encoded image Sfc(r) := s(uj, | r). As a result, 
a typical dMRI data set consists of a collection of such diffusion-encoded images 
{sfc(r)}^j^, whose size K determines the accuracy with which the directions of 
local diffusion flows can be estimated. 

Provided unlimited scanning time, one could measure the diffusion in thou- 
sands of orientations, making it possible to identify the directions of dominant 
diffusion with very high precision. For the reasons explained earlier, however, 
scanning times are always limited, which necessitates restriction of if to a rea- 
sonably small value. This brings us to the central question addressed in the 
present paper: what is a sufhcient number K of diffusion-encoding directions to 
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use? It turns out that, for some realizations of dMRI, the above question can 
be answered in a rigorous manner. In particular, in diffusion tensor imaging 
(DTI) [3}|5|[7{[8| [T6|[T7l , the reconstruction is carried out under the assumption 
that each voxel can support only one diffusion flow as most, which corresponds 
to setting M(r) = 1 for all r in Q. Accordingly, a minimum of if = 7 diffusion- 
encoded images are theoretically sufficient to measure So{r) and recover the six 
non-repetitive components of the symmetric tensors -D(r) by means of least- 
square fitting. (In practice, however, a larger number of gradient directions is 
employed to reduce the estimation variance, with a typical K being between 25 
and 30 18 .) Unfortunately, the accuracy of DTI is known to deteriorate dra- 
matically at the sites where the neural fibres (or bundles thereof) cross, touch 



upon each other, or diverge 13 19 - 27 



The fibre crossing problem in DTI has prompted efforts to develop dMRI 
methodologies which are capable of detecting multiple diffusion flows (or, equiv- 
alently, neural fibre tracts) within a given voxel. One of such techniques is High 
Angular Resolution Diffusion Imaging (HARDI) [T4l[20[|22}|24l[26p8] , which is 



capable of capturing multi-modal diffusion patterns by sampling the spheri- 
cal shell at a much greater number of points (usually between 60 and 100) as 
compared to the case of DTI. Increasing K makes it possible to describe the 
diffusion measurements using much more accurate models. Among these are 
parametric models [29 ■ 34 which allow HARDI signals to be expressed in terms 
of a relatively small number of prototype signal forms. Unfortunately, fitting a 
parametric model often entails minimization of non-convex functionals, which 
is a noise-sensitive and computationally intensive task, prone to the problem of 
local minima. The need to predetermine the optimal number of fitting terms is 
known to be another disadvantage of using the models of the above type. 

The problems associated with parametric modelling of HARDI signals can 
be overcome by using non-parametric models, in which case the signals are 
recovered by projecting the observed data onto properly defined functional sub- 
spaces. In particular, the applicability of spherical Fourier analysis [35] to dMRI 
has been demonstrated in 
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where HARDI signals are approximated by 
truncated series of spherical harmonics (SH). Despite its stability and compu- 
tational efficiency, however, the SH-parameterization involves a relatively large 
number of SHs, which suggests that the SHs cannot be an adequate basis for 
sparse representation of HARDI signals. The main reason for this is rooted in 
the fact that the energies of elementary signals di{u \ r) :— exp {—6 (u-^iDj(r)u)} 
in ([l]) are concentrated alongside the great circles of §^ , whereas the energy of 
SHs is spread all over S^, and, as a result, a relatively large number of SHs are 
needed to effectively "encode" each of di. The inability of the basis of SHs to 
sparsely represent diffusion signals has led to the proposal of spherical ridgelets 
in 



36 37 



where it was shown that it only takes 6 to 8 spherical ridgelets on 
average to represent the HARDI signals with a precision exceeding the precision 
of their representation with 45 SHs. 

The present work takes the ideas of [36 37 one step further and shows that 
the availability of a sparsifying basis for HARDI signals can be used to reduce 
the number of diffusion gradients required for data acquisition. In particular. 
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we suggest to use the theory of compressed sensing (CS) [38 -44 to recover the 
HARDI signals using the number of spherical samples K in a range of values 
typical for DTI (i.e. K e [16, 24]), thus allowing a multi- fibre analysis of dMRI 
data to be performed at the "acquisition cost" of a standard DTI. 

It is worthwhile noting that the ideas of CS have already paved their way into 



the field of diffusion imaging 45 -49 . In this regard, conceptually the closest to 



the proposed approach is the method reported in 46 . In spite of this similarity, 
however, there are two principal distinctions which make the present method 
a more powerful alternative. In particular, the basis functions used in [46] are 
limited to represent an average diffusivity and anisotropy of the white matter, 
thereby neglecting both intra- and inter- voxel variability of tensors Di{r) in 
Q. The ridgelet representation, on the other hand, is a mult iresolut ion tech- 
nique, which possesses an intrinsic ability to deal with a continuum of different 



diffusion scales. Second, the approach in 46 is applied in a "voxel-by-voxel" 
manner and it therefore does not take into consideration the spatial regular- 
ity of diffusion field. The present paper, on the other hand, proposes a novel 
formulation of the problem of CS-based reconstruction of diffusion signals, in 
which the sparsity constraints enforced in the diffusion domain are augmented 
by regularity constraints enforced in the spatial domain. The resulting recon- 
struction problem has the format of a convex minimization problem, which is 



solved using a specially adapted version of the split Bregman algorithm [50 - 52 
As will be shown below, the proposed algorithm results in a particularly advan- 
tageous computational structure which allows the solution to be computed via 
a sequence of simple and easily parallelizable steps. 

The rest of the paper is organized as follows. Section II provides additional 
comments on the input-output structure of the proposed algorithm. The con- 
struction of spherical ridgelets is briefly outlined in Section III, whereas Section 
IV gives a formal description of the proposed reconstruction methodology. Some 
principal details on the numerical implementations of the proposed algorithm are 
summarized in Section V, with the results of our experimental studies reported 
in Section VI. Section VII finalizes the paper with a discussion and conclusions. 



2 Problem Statement 

In the centre of our considerations is the diffusion signal s(u | r) which, when 
normalized by its related 60- image So(r), quantifies the attenuation of MR read- 
out caused by the diffusion of water molecules in the direction u e §^ through 
the spatial position r £ K'^. In practical settings, both u and r are discretized. 
Specifically, restricting u to a discrete set of orientations {u^}^^ prescribes 
the acquisition of diffusion data in the form of K diffusion-encoded images 
{sfc(r)}^^, with each Sfc(r) corresponding to a given u^. In this case, for a fixed 
rg, the vector [si(ro), S2(ro), . . . , Si<-(ro)]"^ G represents a discretization of 
s(u I ro). Note that such a discretization follows a linear measurement model, 
since each sample Sk(jo) can be represented as an inner product of s(u | Tq) with 
a sampling function. In particular, let {<Pfc(u)}^]^ be a Dirac basis of sampling 
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functions defined by 



iPkiu)=Sil-u-Uk), k = l,2,...,K, (3) 

where S denotes the Dirac deha function and the dot stands for the Euchdean 
dot product. Then, formaUy, 

Sfc(ro) = (s(- 1 ro),(/7fc)L2 := / s{u\ro) ipkiu) dr]{u), with /c = 1, 2, . . . , X, 

(4) 

with dr] being the standard rotation invariant measure on S^. 

Next, given a collection of M spherical ridgelets {i^miu)}m=i (defined be- 
low), the signal s(u | r) is assumed to be expandable as 



M 

= ^ c(r)?/'„(u), (5) 



with c(r) e M.^^ being a vector of spherical ridgelet coefficients which depend 
on the spatial coordinate r. It is important to note that the set of spherical 
ridgelets is allowed to be overcomplete, implying dim [Span{-0m(u)}^^]^] < M. 
A practical consequence of this fact is that the definition of coefficients c(r) 
in ([5]) is, in general, not unique. This non- uniqueness is further aggravated 
by the fact that c(r) will have to be recovered from an under-sampled set of 
diffusion measurements, in which case K <C M. Overcoming such a severe 
underdetermination in the problem of estimating the ridgelet coefficients c(r) 
will be possible based on the fundamental premise of the theory of CS, which 
states that an accurate estimation of c(r) is possible if the latter is sufficiently 
sparse and if the sampling and representation bases are sufficiently decorrelated. 



While the sparsity of c(r) is rooted in the very design of spherical ridgelets 37 
the incoherency between the Dirac sampling functions ^ and spherical ridgelets 
stems from the fact that the former have an infinitely small support, whereas 
the latter are "smeared" all over the unit sphere. The above properties of the 
basis of spherical ridgelets yield conditions for an effective application of CS, in 
which case one can obtain a faithful reconstruction of diffusion signals using as 
few as X = 16 diffusion-encoding gradients. 

The proposed algorithm produces an estimate of the ridgelet representation 
coefficients c(r) in ([5]). Once available, the coefficients provide an access to the 
analytical definition of diffusion signals by virtue of ([5]). This can be used in 
a number of ways. One possibility could be to use the ridgelet coefficients to 
compute their associated orientation distribution functions [m] , based on which 
a multi-fibre tractography analysis can be done 53|[54 . Alternatively, ([s]) can be 



used to evaluate the diffusion signals over an arbitrarily fine grid of orientations. 
Subsequently, such refined "measurements" could be fitted using a different 
representation model, whose application to the original data would not have 
been possible without causing severe underestimation errors. Deconvolving the 



refined data to estimate the underlying fibre orientation functions 55 - 58 would 



be another important option to follow. In this paper, we refrain from questioning 
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which of the above possibilities is more advantageous over the others. Our sole 
objective here is to specify a signal processing algorithm which can be used 
to recover HARDI signals, while using the number of diffusion-encoded images 
typical for a standard DTI. 

Finally, it should be noted that the primarily purpose of the proposed 
methodology is to improve the value of HARDI in terms of its time efficiency. 
Since the improvement is achieved through merely decreasing the number of 
diffusion-encoding gradients, the proposed method by no means abrogates the 
use of fast imaging protocols 59 - 6l] to further accelerate the data acquisition. 
Furthermore, an additional speed-up can be achieved via applying CS to recon- 
struct the diffusion encoded images Sfe(r) from their subcritical samples in the 



spectral domain 62 -66 . Generally speaking, we believe this is a combination 
of such software and hardware technologies which will eventually lead to sub- 
stantantial improvements in the practical value of HARDI-based diagnosing. In 
this paper, however, we confine our contribution to showing one particular way 
of attaining this important objective. 



3 Spherical ridgelets 

It is the property of spherical ridgelets to provide sparse representation of dif- 
fusion signals described by ^ which makes them an unparalleled tool for CS- 
reconstruction of HARDI data. To avoid repetitions, in what follows, we present 
only the most principal points of ridgelets design, while their detailed description 



can be found in 37 



Spherical ridgelets are constructed using the fundamental principles of wavelet 
theory 67 68 . Specifically, let x e and p e (0, 1) be a positive scaling pa- 
rameter. Further, let k{x) = exp{—px (x + 1)} be a Gaussian function, which 
we subject to a series of dyadic scalings to result in 

n,{x) = >ii2-^x) = exp [~p J ( J + l) } , (6) 

where j £ N := {0,1,2,...}. Subsequently, the Gaussian- Weierstrass scaling 
function Xj.v : §^ — )■ M at resolution j g N and orientation v € can be 
defined as given by (69 70 



X,,v(u) = ^^^Ac,(n)F„(u.v), VueS^ (7) 

n=0 

where P„ denotes the Legendre polynomial of order n. It is worth noting that 
the L2 energy of Xj.v is concentrated around the spherical point v, with this 
concentration becoming more and more localized when j approaches infinity. 

The spherical ridgelets are designed with the help of the Funk-Radon trans- 
form which, for an arbitrary continuous function / : §^ — )• M, is defined as 

nim^ I /(u)77(u), (8) 
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with i7(v) denoting the great circle perpendicular to direction v, i.e. ct(v) := 
{ugS^|u-v = 0}. Subsequently, following [37|, the semi-discrete frame U of 
spherical ridgelets can be defined as 

U:-{^,-v I ve§2^j = -1,0,1,2,...}, (9) 

where the spherical ridgelet functions i/^j v are obtained from Xj.v according to 

V'j,v(u) = {Xj+i.y - Xj,v} (u), (10) 

with x„i^v(u) = 0. Using ([t]), the ridgelets (10) can be redefined in a closed 
form as (see [35] for details) 

1 °° 2n + 1 

V'j,v(u) = 4^ ('^j+i W - KjH) ■P«(u • v), (11) 

ri=0 

where K_i(n) = 0,V7i and 

f2vr(-l)"/2i:|^, if n is even ^^^^ 
" 1 0, if 71 is odd. 

The set U in (|9]) is infinite-dimensional, and hence is not suitable for practical 
computations. To define a discrete counterpart of U, one has first to restrict 
the values of the resolution index j to a finite set {—1,0, 1, . . . , J}, where J 
defines the highest level of "detectable" signal details. Additionally, the set of 
all possible orientations v G of spherical ridgelets needs to be discretized as 
well. To find a proper discretization scheme, we first note that the construction 



in (111 suggests that the bandwidth of the spherical ridgelets (and therefore the 
dimensionality of the functional space they belong to) increases proportionally 
to 2^ . Since the space of spherical harmonics of degree n has a dimension of 
(n-|-l)2, it seems to be reasonable to define the number of ridgelet orientations 
at resolution j to be equal to Mj = (2^+^mo + 1)^, with mo being the smallest 
spherical order resulting in Ko(mo) < e for some predefined < e <C 1 (e.g. 
e = 10~6). Consequently, for each j, a total of Mj orientations {v^}f£^^ are 
chosen so that a discrete counterpart of U can now be defined as 

Vd = [>P,^^. Ij = -1,0, 1, . . . J, z = 1,2, . . . ,M,} . (13) 

where the subscript d stands for "discrete". It should be noted that, although 
the set Urf is composed of continuously defined functions, its dimension is finite, 
since Ud consists of a total of M = X]j=-i(2"'^^™o + 1)^ spherical ridgelets. To 
slightly simplify our notation, in what follows, the spherical ridgelets in Vd will 
be indexed as ?Am(u), with m = 1, 2, . . . , M being a combined index accounting 
for both different resolutions and orientations. 

Given a sampling set of K diffusion-encoding orientations {ufcj^j^, one can 



use (111 to compute the values of the spherical ridgelets in over the sampling 
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sei[2 The resulting values can be stored into a. K x M matrix A defined as 

-01 (Ul) V'2(ui) ... VAf(ui) 
V'1(U2) V'2(U2) ... V'Af(u2) 

_ i/'i(uk) V2(uk) ... Vm(uk) _ 

Then, for a given vector s(r) :— [s(ui | r), s(u2 | r), . . . , s{uk \ G of the 
measured values of a diffusion signal s(u | r) at the spatial location r, the model 
([5]) asserts the existence of representation coefficients c(r) € M.^ such that 

s(r) =v4c(r) + e(r), (15) 

where e(r) accounts for both model and measurement noises. Clearly, the non- 
negligibility of e(r) along with the fact that K <^ M makes the problem of 
recovering the representation coefficients c(r) from s(r) a very challenging in- 
verse problem, our solution to which is presented next. 



4 Proposed reconstruction framework 

Let n represent the volume within which diffusion measurements are acquired. 
Also known as a region of interest, is assumed to be a bounded rectangular 
subdomain of K^, i.e. := [0, L^] x [0, Ly] x [0, L^] C K^. Let fid be a discrete 
subset of ri, which represents the spatial locations at which the diffusion signal 
is measured. Specifically, Qd is assumed to be a uniform lattice which can be 
formally defined as 



— {xii,yi2, Zi^} ^ii — ^L^, ~ —Ly, ~ —Lz\ , (16) 



where < ii < N^, < 12 < Ny, and < 13 < iV^ are sampling indices in the 
direction oi x, y and z coordinates, respectively. 

Let K be the number of diffusion-encoding gradients used for HARDI data 
acquisition, and let the corresponding gradient orientations be denoted by {u^}^ 
where u^. S S^. For each of these values of u^., MRI measurements result in 
its associated diffusion-encoded image Sk, which can be formally viewed as a 
mapping from Sl^ to M. For the sake of convenience, each Sk can be stored 
and manipulated as an x Ny x Nz array of real numbers, namely Sk € 
^N^xNyxN^ _ Alternatively, at a given coordinate r G fid, one can combine the 
values Si(r), S2(r), . . . , s/^ (r) into a column vector s(r) := [si(r), 52(1"), ... , s^f (r) 



(as it was already done in ( 15 )). This vector can then be regarded as a vec- 
tor of discrete measurements of an associated HARDI signal s(u | r) measured at 
orientations {u^j^j^. It is worth noting that, according to the above notations. 



^ Since the definition in involves an infinite summation, the latter needs to be truncated 
to render the computations practical. In practice, we truncate the summation to index Umax 
for which the magnitude of the summand drops below 10~®. 
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the value Sk (r) admits a twofold interpretation, viz. either as the fc*^ coordinate 
of vector s(r) or as the value of image Sk at spatial position r. 

When combined together, the continuum of vectors s(r) can be regarded as a 
discrete vector field s : M^, in which case s(r) has a natural interpretation 
of the value of s corresponding to position r. The vector space 9J of such vector 
fields can be endowed with the standard inner product 



{s\s'U = E W^-^(r) = E E sUr) = E(4, 4), (17) 

rGfid r^Qd k—1 



fe=l 



with {si, si) = J2ren 'Sfc(r) being the standard inner product between the 



scalar- valued images s], and s^. Accordingly, congruent to the definition in (17) 
the £2-noTm of s e 23 is defined as 



[ E 



1/2 



K 



fc=l 



1/2 



(18) 



where || • II2 and || • \\p denote the Euclidean vector and Frobenius matrix norms, 
respectively. 

Another norm on 2J that we shall make use of is the total variation (TV) 
semi-norm which is defined as follows. First, let us define the total variation of 
the fc"^ component Sk of the field s in a standard manner as 



|Sfc||TV 



E 

rend 



pec(r) 



1/2 



(19) 



where C (r = (x;, , , z^g)) = {{x^^^i,yi^, z^^), {x,^,yi^^i, z,^), {xi^.y.^^, z^^^i)} 
is a 3-neighbourhood (causal) clique of voxel r. Consequently, the TV norm 
of the discrete vector field s can be defined in terms of the TV-norms of its K 
components as 

K 

Plb,TV= [Ell^fcllTv] (20) 

k = \ 



Thus, for example, a — 2 was used in the TV-denoising method reported in 71 
In this paper, we use a = \. 



Now, let {'0m}m=i * set of spherical ridgelets defined by (13), which is 
assumed to be rich enough so that each HARDI signal can be expressed accord- 
ing to ([5|. Analogously to the discrete measurements s(r), the representation 
coefficients c(r) corresponding to different voxels r can be aggregated into a 
vector field c e il, where 11 : fid — (with c(r) being the value of c observed 
at r). Although it is possible to endow the vector space il with both the £2- 



and TV- norms by analogy with (18) and (20), it will be particularly useful to 
consider the ^i-norm of c which can be defined as 



iciiu.i- E ii^wiii= EEic'^wi- 



(21) 



r&ld fe=l 
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Using the definitions of tfie vector fields QJ and il as well as the definition 



of A in (14), a connection between QJ and il is established by means of a linear 



map ^ : il -> 93 that is defined as given by 

: il^ QJ : c(r) ^ s(r) = Ac{r), Vr e ^a. (22) 

Consequently, using A, one can define the HARDI data formation model as 

s ^ A{c} + e, (23) 

where e € QJ is supposed to account for both measurement noise and modelling 
errors, and it is assumed to have a relatively small €2-uorm ||e||2j^2 < £• 



The model ( 23 ) suggests a reduction of the problem of estimation of HARDI 
signals to the problem of estimation of their corresponding representation coeffi- 
cients c from the discrete and noisy measurements s. Furthermore, as our main 
intension has been to recover the coefficients c using as few diffusion-encoding 
gradients as possible (implying K ^ M), there is an infinite number of solu- 
tions which would fit the constraint ||.4{c} — s||(jj^2 < £■ However, if it is known 
a priori that, for each r e fid, the vector of representation coefficients c(r) is 
sparse, then a useful estimate of c can be obtained as a solution to the following 



convex optimization problem 38 - 44 



minllcllai.i (24) 

C 

s.t. P{c}-s||qj,2 < e. (25) 



It should be noted that the optimization problem ( 24 ) is equivalent to solving 



mm 2^ ||c(r)||i (26) 

s.t. J2 \\Ac{r)-.s{r)\\l<e^, (27) 
rend 

and therefore, under the assumption of spatially homogeneous noise e, the prob- 



lem (24) is separable in the spatial coordinate r. This means that an optimal 



field c can be recovered by solving for its components 

min||c(r)||i (28) 

c(r) 

S.t. M{c(r)} - s(r)||2 < (N^NyN^y'/h, (29) 

independently at each r e 57^. 

While computationally attractive, the above solution is suboptimal, since it 
completely disregards the dependencies which are likely to exist between spa- 
tially adjacent HARDI signals. A possible way to take such dependencies into 
consideration is to require the noise-free version of the measured signal s to 



possess a minimal TV norm among all possible candidate solutions 72 . This 
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requirement can be translated into the following minimization problem 



mm |||c||2jj+7 ||yl{c}||-iT,Tv} 

s.t. \\A{c} - s||,jj,2 < e. (30) 
where the role of 7 > is to balance the relative influence of the sparse and 



TV terms in the above cost function. The optimization problem ( 30 1 can be 
rewritten in its equivalent Lagrangian form 



(31) 



for some optimal values of A > and ^ > 73 . 

Below, we are going to specify a particular, computationally efficient method 



instances of (31 1 



for solving (31 1. In this connection, it is instructive to outline the following two 



4.1 Sparse-only reconstruction 



When /X = 0, the functional in (31 1 becomes separable in the spatial variable r 



in the sense that, in such a case, an optimal c can be recovered by solving 



mm 

c(r) 



1, 



Ac{r) - s(r) 



A||c(r) 



(32) 



for each c(r) independently. Note that ( 32 ) can be considered to be a Lagrangian 
form of the optimization problem ( 28 1 . There exist a broad spectrum of methods 



which could be used for solution of (32 ). Some particularly attractive algorithms 



seem to be those exploiting the principle of iterative shrinkage (aka iterated 
thresholding) 74-77 . While the non-differentiability of £i-norm in (|32|) rules 



out the applicability of gradient-based optimization tools, iterative shrinkage 



methods are capable of finding a solution of ( 32 ) through iterative application 
of a first-order, fixed-point update rule. Specifically, many algorithms of this 
kind find an optimal solution as a stationary point of the sequence of estimates 
produced by 



c(r) 



t+i _ 



5a/. {c(r)* + (s(r)-Ac(r)*)}, 



(33) 



where Sr{t} = sign{t){\t\ 
V is chosen to obey v > 



— r)_|_ denotes the operator of soft thresholding and 
In the present paper, a modification of the 
iterative update in (33), known as the jast iterative shrinkage thresholding al- 
gorithm (FISTA) f78f7 was employed due to its considerably faster convergence 
as compared with many alternative "accelerated" methods. 

It should be emphasized that, while being suboptimal from the viewpoint 



of spatial-domain regularity, the solution of ( 32 ) through iterative shrinkage is 



advantageous in two important practical ways. First, it suggests considerable 



12 



storage reduction, since the thresholding operator in ( 33 1 sets to zero the rep- 
resentation coefficients with amplitudes less or equal to \jv in absolute value. 
It makes it possible to use sparse data formats to store and manipulate the 
representation coefficients. Second, the fact that the estimation of c(r) is per- 
formed at each voxel independently suggests a natural way to speed up the 
overall estimation process though parallel computing on a multicore system. 

4.2 TV-only reconstruction 



When A = 0, solving the optimization problem (31) is equivalent to simultane 



ously solving K optimization problems of the form 

^|IWc}],-.fc|l|. + M||[^{c}]fc||Tv}, (34) 

where fc = 0, 1, . . . , if — 1 and [y^{c}]fc denotes the fc-th component of the vector 
field A{c\ G QJ. Let [yl{c}]fc be denoted by u^, i.e. ■— [A{c}]k- Then, 



reformulated with respect to Ufc, the problem (34) can be rewritten as 

1 



mm 



^u'^k - SkWp + l-J-WukWrv } , (35) 



in which case it can be recognized as the problem of TV-denoising of the 



diffusion-encoded image st 72 . It is important to note that, in contrast to 



( 34 ) , the problem ( 35 ) can be solved for each k independently, in which case we 



say that the estimation becomes separable in the diffusion direction. 



The current arsenal of methods which can be used for solving ( 35 ) is broad. 
Originated from the work of Rudin et a/ 1 72 , TV-denoising methods currently 
include gradient-based optimization methods [79 , first-order methods |80j, it er- 
ative shrinkage methods 81 82 , and Bregman-type iterative algorithms [50[[5T| 
[831 . The first-order methods appear to be a particularly attractive option when 
relatively large data arrays need to be processed, which is the case relevant to 



dMRl. Consequently, in the present paper, we employ the algorithm of 84 for 
the simplicity and elegancy of its implementation as well as for its outstanding 
convergence properties. 

5 Solution using Split Bregman Algorithm 

Directly solving the original problem (IsTl) is difhcult because of the compound 



nature of the regularization it involves. The split Bregman approach 50,51 



allows reducing (31) to a simpler form through introduction on an auxiliary 



variable u € 03, which can be viewed as a noise-free version of the data field s. 



Particularly, using u (31 1 can be redefined as 

min|^||M-s||^2+-^l|c||H.i + ^J■\\u\\<a,Tv\ (36) 

c,u L Z J 

s.t. M{c} - m|||j,2 = 0. 



13 



Then, starting from an arbitrary 6° € QJ, the Bregman algorithm 85 finds 
optimal c and u through the following iterations 



(||2 
23, 



c*+i) = arginin - s\\% .,+\ ||c||u,i + \\u\\y^^^ + i"" " "^^^^ " ^'H 

(37) 

for some 7 > The functional in ( |37| is supposed to be minimized over 
two variables, i.e. u and c. However, due to the way the ii and TV compo- 
nents of this functional have been split, the minimization can now be performed 
efficiently by iteratively minimizing with respect to u and c separately. The 
resulting iteration steps are 

Stepl: c*+i =argmm{^P{c}-(u*-6*)||2,2 + A||c||ua} (38) 

Step 2: = argmm{i||M - s|||, 2 + ^H^* - (-4{c*+^} + b')\\l,2 + I|w||w,tv}. 

Note that the functional at Step 2 contains two quadratic terms which can be 
combined together to result in 

o t+l • rl+7ii s + 7(yl{c*+i} + 5*)||2 „ „ 1 

step 2: u ^ =argmin| ^ \\u — jj^^j ^ + M ll'"ll'»,Tv|- 

^ (39) 



To cause a substantial reduction in the value of the cost functional in (37), 
Step 1 and Step 2 should be applied recursively for a predefined number of 
times before the Bregman parameter 5* is updated according to (37). It was ar- 



gued in 51 , however, that the extra precision gained through such a repetitive 
application of Step 1 and Step 2 is likely to be "wasted" when 6* is updated. 



Consequently, it was suggested in 51 to perform these steps only once per itera- 
tion cycle. It is interesting to note that, in this case, the split Bregman algorithm 
transforms into the alternating directions method of multipliers (ADMM) |52| , 
whose convergence is guaranteed by the Eckstein-Bertsekas theorem [86] (see 
also Theorem 3.1 in [52^). 

The final algorithm is summarized below. Lines 3-4 of Algorithm 1 corre- 
spond to Step 1 in ( [38| , while lines 5-6 correspond to Step 2. An even more 
important fact to notice is that the optimization problem in line 4 is separable 
in the spatial coordinate r. This optimization, therefore, can be performed at 
each voxel independently as discussed in Sect ion [4T| Moreover, the optimization 
problem in line 6 is separable in the diffusion coordinate fc, and hence it amounts 
to applying TV-denoising to each of the K components of u independently as 
discussed in Section 1421 



■^Note that the algorithm is guaranteed to converge for any 7 > 0. In this work we use 
7 = 0.5. 
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Algorithm 1 ADMM algorithm for sparse- TV reconstruction of HARDI signals 

1: b ^ 0, u ^ s 

2: while "c keeps changing" do 
3: d ^ u — b 

4: c argminc llll^jc} - d\\% .^ + ^ l|c||ii.i} 

5: d-^ (l + 7)-i(s + 7(yt{c} + 6)) 

6: argminu jjHu - d|||,^2 + ll"ll'»,Tv| 

7: 6-^6+ {A{c} - u) 

8: end while 
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Figure 1: Phantom #1: (Upper row of subplots) The orientations of the indi- 
vidual diffusion flows and their combination; (Lower row of subplots) Examples 
of the resulting (noise-free) diffusion-encoding images corresponding to four dif- 
ferent diffusion-encoding directions. 
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Figure 2: Phantom #2: (Upper row of subplots) The orientations of the indi- 
vidual diffusion flows and their combination; (Lower row of subplots) Examples 
of the resulting (noise-free) diffusion-encoding images corresponding to four dif- 
ferent diffusion-encoding directions. 



SNn = i SNR = 24dB SNR = 18dB SNH=12dB 




Figure 3: (Upper row of subplots) Diffusion-encoding images of Phantom #1 
corresponding to u = [1, 1, 1]/a/3 and SNR — oo, 24, 18 and 12 dB; (Lower 
row of subplots) Diffusion-encoding images of Phantom ^2 corresponding to 
the same u and the same values of SNR. 
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6 Results 



6.1 Technical details of the experimental study 

Both the choice of difFusion-encoding directions {u.k}^^i and of the orientations 
of spherical ridgelets require the use of a samphng scheme. In dMRI, one of the 
standard methods used to distribute a given number of spherical points over 
in a quasi-uniform manner is by means of the electrostatic repulsion algorithm 
(also known as the Thomson problem). However, since the diffusion signals are 
symmetric (implying s(u | r) = s(— u | r)), it is a unit hemisphere, not the entire 
S^, which actually needs to be discretized. In view of the absence of a formu- 
lation of the Thomson problem for hemispherical domains, a common practice 
is to run the standard procedure for twice as many points as needed, followed 
by keeping only a half of the resulting configuration. However, as the retained 
points are not explicitly constrained to lie on a hemisphere, they may include 
nearly antipodal pairs which are likely to introduce undesirable dependencies 
between the diffusion measurements as well as between the basis functions. This 
limitation can be overcome through adapting a different sampling strategy. Par- 
ticularly, in this paper, both the diffusion-encoding directions and the orienta- 
tions of spherical ridgelets have been defined by using the method of generalized 



spiral points 87 88 , in which the sampling points are arranged along a spiral 
in such a way that the distance between the points along the spiral is approxi- 
mately equal to the distance between its coils. This method is easily adaptable 
for sampling of the "northern" hemisphere (i.e. {u S | u • [0, 0, 1]"^ > O}), 
providing a nearly uniform, unique and analytically computable coverage which 
is in no respect inferior to the one produced by solving the Thomson problem. 

To assess the performance of the proposed algorithm under controllable con- 
ditions, experiments with simulated data sets have been performed. In this case, 
the HARDI signals were generated according to model ([ij with different values 
of M{r), Di{r), and so(r) = 1, Vr. The resulting signals were contaminated 
by variable levels of Rician noise, giving rise to a set of different SNRs. In this 
work, we adapt the standard definition of the SNR as 

SNR = 201ogiof^nr^)' (40) 
where s and s denote an original signal and its noise-contaminated version, re- 



spectively, and the norms are computed as defined by (18). It is worthwhile 
noting that the optimal values of regularization parameters A and fj, in (31 1 
are normally a function of the noise level. In the present paper, however, no 
attempts have been extended to optimize these values for different SNRs. In- 
stead, it was found that A = 0.03 and ^ = 0.05 provided acceptable estimation 
results in all the simulation scenarios, and hence these values have been used 
throughout the whole study. 



Following J37 



, the scaling parameter p in (|6| was set to 0.5 and the resolution 



parameter J in (13) was set to be equal to 1, corresponding to a total of three 



resolution levels. The number of spherical ridgelet orientations were predefined 
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with niQ — 4, resulting in M_i = 16, Mq — 49 and Mi = 169 ridgelets spanning 
the resolution levels j — —1, j — and j = 1, respectively. Thus, the total 
number of spherical ridgelets used in the reconstruction was equal to 234. 

To quantitatively compare the reconstruction results produced by the pro- 
posed and references methods for different numbers of sampling directions K 
and various SNRs, three performance measures were used. The first of the three 
was the normalized mean-squared error (NMSE) defined as 

with s(r) being a reference HARDI signal corresponding to location r and s(r) 
being its estimate. Depending on the nature of a specific experiment, the refer- 
ence signal can be either a simulated signal discretized at 642 spherical points 
obtained by the 3rd order tessellation of the icosahedron or a signal recon- 
structed using a maximum possible number of diffusion-encoding orientations. 

One of the most valuable outcomes of HARDI is in providing an access to 
computation of orientation distribution functions (ODFs) - the functions whose 



modes are likely to coincide with the direction of local diffusion flows 14 . Both 



the SH-based |28| and ridgelet-based 37 methods of reconstruction of HARDI 



signals come with analytical expressions which relate the HARDI signals to their 
corresponding ODFs. The latter can in turn be used to recover the directions 
of local diffusion flows (or, equivalently, the orientations of their related fibre 



tracts) using, e.g., the steepest ascent procedure detailed in 37 . Suppose Uq is 
the true direction of a diffusion flow and u is its estimate. Then, the angular 
orientation error S can be defined (in degrees) as 

= arccos(uo • u). (42) 

TT 

In this paper, as a performance measure, we use an average angular orientation 
error which is obtained by averaging the values of S computed for all "fibres" 
within a specified Qd- 

The last performance measure used in this work is the probability Pd of false 
fibre detection. To define Pd, let M{r) be the true number of fibre tracts passing 
through voxel r (as defined by model ([ij). Also, let M{r) be an estimated 
number of fibres, which is equal to the number of modes (maxima) of the ODF 
recovered at position r. Then, one can define 



1 ^ |Af(r)-Af(r)| 



N^NyN. rt^, Mir) 



■ 100%. (43) 



In addition to the quantitative comparison, the reconstruction results will 
be evaluated through visual comparison as well. In this paper, we choose to 
visualize spherical functions by means of 3-D surface plots. Such a plot tends 
to project away from the origin of in the directions along which a spherical 
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function is maximized, while passing near the origin in the directions where the 
function approaches zero. 

Finally, our choice of reference methods was motivated by the scope of the 
main statements made in this paper. First, since we argue that the frame of 
spherical ridgelets is optimally suited for CS-based reconstruction of HARDI sig- 
nals, its performance has to be compared with that of alternative representation 
systems. In particular, the basis of spherical harmonics up to the order 8 inclu- 
sive has been used for a different definition of the sensing matrix A in ( 14 ) . (Note 
that, in the case of a real and symmetric analysis, this SH-basis consists of 45 
functions.) Additionally, the representation system proposed in 46p7| has been 
exploited in the comparative study as well. This system is formed by applying a 
set of rotations to a Gaussian kernel of the form d{u) = exp{— 6 {u^Dq u)}, with 
b defined as in ([I]) and Dq equal to a (diagonal) diffusion tensor having a mean 



diffusivity of 766 mm^/s and a fractional anisotropy of 0.8. Following 47 , the 
number of rotations (and hence the number of Gaussian basis functions) was set 
to be equal to 253. For the convenience of referencing, the CS-based reconstruc- 
tion methods using the spherical ridgelets, the 8th-order spherical harmonics, 
and the rotated Gaussian kernels will be referred below to as the RDG, SH8 
and GSS algorithms, respectively. 

To assess the significance of the proposed spatial regularization, all the above 
algorithms have been applied with two different values of /i in (31 ), viz. /i = 
and /i = 0.05. Note that, in the first of these cases, the spatial regularity is 
ignored, which leads to the sparse-only reconstruction discussed in Section [4T| 
In the second case, on the other hand, the spatial regularity is taken into account 
and the reconstruction is performed by means of the split Bregman algorithm 
of Section [5] 



6.2 In silico experiments 

To assess the performance of the proposed and reference methods under con- 
trollable conditions, two simulated data sets were used. The first set (referred 
to below as Phantom #1) had a spatial dimension of 12 x 12 pixels, and con- 
sisted of two "fibres" crossing each other at the right angle as it is shown in 
the upper row of subplots of Fig. [T] In addition, each pixel in the set was as- 
signed an extra diffusion flow in the direction perpendicular to the image plane. 
As a result, the number of diffusion components A/(r) in Phantom #1 varied 
between 1 and 3. Subsequently, model ([!]) was used to generate corresponding 
diffusion-encode images {sk}^^i for a range of different values of K. Two dif- 
ferent values of 6, namely b — 1000 s/mni^ and b = 3000 s/mm^ were used for 
data generation. The diffusion tensors Di{r) in ([l]) were obtained by applying 
rotations to a tensor of the form Dq — diag ([a, /3, /3]), where a and /3 were equal 
to 1700 • 10"^ and 300 • 10"^, respectively. Note that the mean diffusivity and 
fractional anisotropy of Dq are equal to 766 mm^/s and 0.8, respectively. Thus, 
the same diffusion tensors were used for data synthesis and for the construction 
of basis functions in the GSS algorithm, thereby allowing the latter to perform 
under the best possible conditions. 
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Figure 4: (Upper subplot) Original ODFs of Phantom #1; (Middle row of sub- 
plots) The ODFs recovered by the SH8-CS, GSS-CS, and RDG-CS algorithms, 
respectively; (Bottom row of subplots) The ODFs recovered by the SH8-TV, 
GSS-TV, and RDG-TV algorithms, respectively. 
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Figure 5: (Upper subplot) Original ODFs of Phantom #2; (Middle row of sub- 
plots) The ODFs recovered by the SH8-CS, GSS-CS, and RDG-CS algorithms, 
respectively; (Bottom row of subplots) The ODFs recovered by the SH8-TV, 
GSS-TV, and RDG-TV algorithms, respectively. 
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The lower row of subplots in Fig. [T] depict four examples of the diffusion- 
encoding images obtained for Phantom ^1 before their contamination by Rician 
noise. One can see that the images are piecewise constant functions, which 
appears to be in a good agreement with the bounded- variation model suggested 
by (31 ). However, real images may be more complicated than that. Accordingly, 
to test the robustness of the proposed regularization scheme, a different in silico 
phantom was designed. Phantom #2 had a spatial dimension of 16 x 16 pixels 
and it was obtained through supplementing the configuration of Phantom #1 
by an additional circular "fibre" as shown in the upper row of subplots in Fig. [2j 
The lower row of subplots of the figure show a subset of the resulting diffusion- 
encoded images, which can be seen to no longer exhibit a piecewise constant 
behaviour characteristic for Phantom #1. 

The simulated diffusion-encoded images were contaminated by three differ- 
ent levels of Rician noise, giving rise to SNR of 24, 18 and 12 dB. Some typical 
examples of the resulting images are demonstrated in Fig. [3] where the up- 
per row of subplots depict a noise-free version of one of the diffusion-encoded 
images pertaining to Phantom ^1 along with its noise-contaminated counter- 
parts. The lower row of subplots in Fig. |3] depict an analogous set of examples 
for Phantom #2. Observing the figure, one can see that the SNR values have 
been chosen so as to cover a range of possible noise scenarios, which could be 
characterized as moderate-to-severe contamination. 

As it was mentioned earlier, in our in silico study we compared the perfor- 
mance of three different representation bases, i.e. spherical harmonics (SH8), 
Gaussian kernels (GSS) and spherical ridgelets (RDG). All the resulting al- 
gorithms have been further subdivided into two different types, depending on 
whether or not the spatial regularization was engaged. Thus, in the absence of 
the spatial regularization (corresponding to ^ — 0), the reconstruction has been 



performed on a voxel-by- voxel basis, as detailed in Section 4.1 For the conve- 



nience of referencing, the corresponding algorithms will be referred to below as 
SH8-CS, GSS-CS, and RGD-CS. In the case of ^ > 0, the estimation has been 
carried out using the split Bregman method of Section [5] The corresponding 
algorithms will be referred below as SH8-TV, GSS-TV, and RGD-TV. 

The upper subplot of Fig. [4] shows the original field of ODFs of Phantom #1 
(corresponding to 6 = 3000 s/mm^), which have been computed based on Tuch's 
approximation [14| (i.e. by applying the Funk-Radon transform to the diffusion 
signals). At the same time, the middle row of subplots of Fig. [4] show the ODFs 
recovered by (from left to right) SH8-CS, GSS-CS and RDG-CS with K = 16 
and SNR = 18 dB. One can see that the inability of the SH basis to sparsely 
represent HARDI signals results in a poor performance of SH8-CS. A better 
result is obtained with GSS-CS, which uses a basis of rotated Gaussian kernels, 
and therefore has a potential to represent the HARDI signals in a sparse manner. 
Unfortunately, the excessive correlation between the Gaussian basis functions 
adversely affects the ability of this method to withstand the effect of noise. 
Consequently, the reconstruction obtained using GSS-CV suffers from sizeable 
errors. The RDG-CS method, on the other hand, provides an estimation result 
of a much higher quality, albeit some inaccuracies are still noticeable in the 
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Figure 6: NMSE obtained using the compared methods for different phantoms, 
SNRs and 6- values. 
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Figure 7: Average angular error 6 obtained using the compared methods for 
different phantoms, SNRs and 6-vahies. 
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Figure 8: The rate of false fibre detection Pd obtained using the compared 
methods for different phantoms, SNRs and 6-values. 
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central part of the phantom. The reconstruction accuracy improves dramatically 
when the spatial regularization is "switched on" , as it is demonstrated by the 
bottom row of subplots in Fig. |4j Specifically, while SH8-TV is still unable to 
provide a valuable reconstruction, the estimates obtained using GSS-TV and 
RDG-TV represent correctly the "flow structure" of Phantom #1. Moreover, 
among the latter two methods, RDG-TV is clearly the best performer, resulting 
in a close-to-ideal recovery of the original ODFs. The superiority of RDG- 
TV over the alternative methods is further evident in the results presented by 
Fig. [5] which depicts the reconstructions obtained for Phantom #2 (with the 
same values of b, K and SNR as above). 

In general, the reconstruction results obtained using SH8-CS and SH8-TV 
have been observed to be of a lower quality in comparison to the other methods 
under consideration. For this reason, in what follows, only the GSS and RDG 
methods are compared. Thus, Fig. |6] contrasts the performances of GSS-CS, 
GSS-TV, RDG-CS and RDG-TV in terms of the NMSE criterion. One can see 
that the best performance here is attained by the RDG-TV algorithm, which 
results in the smallest values of NMSE for both phantoms and for all the tested 
values of 6, SNR and K. It is also interesting to note that the incorporation of 
spatial regularization allows GSS-TV to outperform RDG-CS, with the effect 
of the regularization becoming more pronounced at lower SNRs. On the whole, 
all the NMSE curves demonstrate an expected behaviour, with the error values 
increasing proportionally with a decrease in SNR, while going down with an 
increase in the number of diffusion-encoding gradients K. However, as opposed 
to the others, the NMSE curves obtained with RDG-TV are characterized by 
a relatively low rate of convergence, which indicates a reduced sensitivity of 
RDG-TV to the value of K. 

The above algorithms have been also compared in terms of the angular error 
(42). The results of this comparison are summarized in Fig. [7j which again 
indicates that the most accurate reconstruction is obtained using the RDG-TV 
method. In general, the angular error tends to grow as SNR decreases and 
to converge to a minimum as K increases. As opposed to the case of NMSE, 
however, there is an additional dependency of the angular error on the type of 
a phantom in use as well as on the &-value. In particular, the errors obtained 
for Phantom #2 are (on average) greater than those obtained for Phantom #1. 
This discrepancy is rooted in the fact that Phantom #2 has a more complex 
"fibre structure" as compared to Phantom ^1. In particular, while the "fibers" 
of Phantom #1 are designed to cross each other at the right angle, the "fibres" of 
Phantom #2 are allowed to decussate at much smaller angles, which makes them 
much harder to resolve. Moreover, this effect becomes more noticeable with a 
decrease in the 6-value, which reduces the resolution of q-ball imaging. Finally, 
we notice that, on average, GSS-TV performs better than RDG-CS (though still 
worse than RDG-TV), which justifies the value of spatial regularization. 

The comparison in terms of the rate of false fibre detection Pd ( 43 1 was last 
in the line of our in silico performance tests; its results are shown in Fig. [Sj 
One can see that, in the case of Phantom ^1, RDG-TV yields a virtually zero 
false detection rate for both values of b, whereas the other methods result in 
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considerably higher values of Pd (mainly due to the detection of spurious local 
maxima in the estimated ODFs). The situation is different for Phantom #2, 
where all the compared methods yield sizeable errors (especially for h = 1000 
s/mm^). However, in comparative terms, the most accurate reconstruction is 
still obtained by means of the proposed RDG-TV algorithm. 



6.3 In vivo results 

As the next validation step, experiments with real HARDI data were carried 
out. The proposed algorithm was tested on human brain scans acquired on 
a 3-Tesla GE system using an echo planar imaging (EPI) diffusion-weighted 
image sequence. A double echo option was used to suppress eddy-current re- 
lated distortions. To improve the spatial resolution of EPI, an eight channel 
coil was used to perform parallel imaging by means of the ASSET technique 
with a speed-up factor of 2. The data were acquired using 51 gradient direc- 
tions (quasi-uniformly distributed over the northern hemosphere) with b = 1000 
s/mm^. In addition, eight baseline (60) scans were acquired, averaged and used 
for normalization. The following scanning parameters were used: TR — 17000 
ms, TE = 78 ms, FOV — 24 cm, 144 x 144 encoding steps, and 1.7 mm shce 
thickness. All scans had 85 axial slices parallel to the AC-PC line covering the 
whole brain. 

The main question addressed through the in vivo experiments has been 
whether or not it is possible to supersede the spatial regularization by pre- 
filtering of HARDI signals. To this end, the RDG-CS algorithm was applied 
first to the HARDI data containing the full set of X = 51 diffusion gradients. 



(Note that such dense reconstruction is analogous to the one reported in 37 



where the latter is shown to outperform the SH-based estimation 28 .) The re- 
sulting ODFs have been used as a fiducial against which different reconstruction 
results were compared. 

As the next step, three different subsets of 16, 24 and 32 spherical points 
were composed out of the original set of 51 diffusion gradients. Within each of 
these subsets, their corresponding points were chosen so as to result in a quasi- 
uniform coverage of the northern hemisphere. Accordingly, the HARDI data 
were rearranged into three data sets of size 144 x 144 x 85 x 16, 144 x 144 x 
85 X 24 and 144 x 144 x 85 x 32 to emulate compressed sensing data acquisition. 
The above sets were used to assess the performance of different reconstruction 
methods. Unfortunately, we have not succeeded to find conditions under which 
the SH8 and GSS algorithms would provide stable reconstruction results (either 
with or without pre-filtering) . For this reason, only the RDG-CS and RDG-TV 
algorithms are compared below. 

The upper row of subplots in Fig.|9]show the generalized anisotropy (GA) [l4j 
image of a coronal cross-section of the brain along with the reference field of 
ODFs corresponding to the region indicated by the yellow rectangular. Anatom- 
ically, this region is expected to contain the fibre bundles of corona radiata as 
well as those of superior longitudinal and arcuate fasciculi. The middle row of 
subplots in the same figure depict the ODFs reconstructed by RDG-CS using 
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Figure 9: ( Upper row of subplots) A coronal GA image and the ODF field 
of the indicated region recovered by RDG-CS with K = 51; {Middle row of 
subplots) Estimated ODF fields obtained using RDG-CS with X = 16, X = 24 
and K = Z2; [Bottom row of subplots) Estimated ODF fields obtained using 
RDG-TV with K = 16, K ^ 24 and K = 32. 
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Figure 10: ( Upper row of subplots) An axial GA image and the ODF field 
of the indicated region recovered by RDG-CS with K = 51; {Middle row of 
subplots) Estimated ODF fields obtained using RDG-CS with K = 16, K = 24 
and K = 52; {Bottom row of subplots) Estimated ODF fields obtained using 
RDG-TV with K = 16, K ^ 24 and K = 32. 
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Table 1: NMSE computed between the dense and CS-based reconstructions 
o btained with RDG-CS and RDG-TV 





Pertaining to Fig. 9 


Pertaining to Fig. 10 


K = 16 K = 2A K = 32 


K =16 K = 24 K = 32 


RDG-CS 
RDG-TV 


0.097 0.064 0.043 
0.022 0.011 0.003 


0.091 0.053 0.037 
0.018 0.009 0.002 



K = 16, 24 and 32 diffusion gradients. One can see that the quahty of recon- 
struction progressively improves as K increases. It is important to note that, 
before applying the RDG-CS algorithm, the diffusion-encoded images had been 
pre-processed by a TV filter to reduce the effect of measurement noises on the 
estimation result. However, this pre-processing appears to be not nearly as 
effective as the spatial regularization of the RDG-TV algorithm, whose recon- 
struction results are shown in the bottom row of subplots in Fig. |9] The above 
conclusion is further supported by an additional example of Fig. |10[ which shows 
the reconstructions pertaining to the indicated area within an axial cross-section 
of the brain. (The relevant fibre bundles here are those of cingulum and corpus 
callosum). As in the previous example, one can see that the most accurate 
reconstruction is attained by means of the proposed RDG-TV method. The 
superiority of RDG-TV is also confirmed by the quantitative figures of Table [T] 
which summarizes the NSME obtained by the compared algorithms for different 
values of K. 

7 Discussion and Conclusions 

When considered as a whole, the HARDI signals which pertain to a given volume 
of interest can be described as multi- valued (or, more generally, measure- valued) 
fimctions from a subset of M.^ to the space of square-integrable spherical func- 
tions L2(§^). Such functions can be thought of as if they had two "modes of 
variation" - one in the spatial and another in the diffusion domain. Although 
applying various inverse problems (a particular instance of which is addressed 
by the theory of CS) along the spatial and diffusion coordinates independently 
is not new to the community of medical imaging scientists, formulating a CS 
reconstruction problem in both domains simultaneously has not been proposed 
before. Accordingly, the present paper introduced the RDG-TV algorithm which 
exploits the above idea and can be used for reliable reconstruction of HARDI sig- 
nals from as few as if = 16 diffusion-encoded scans (as compared to 60-80 scans 
required by existing reconstruction tools). The algorithm exploits the fact that 
HARDI signals can be sparsely represented by spherical ridgelets in the diffusion 
domain, while their associated diffusion-encoded images have bounded variation 
in the spatial domain. Moreover, it has been shown experimentally that either 
using different representation bases or excluding the spatial regularization would 
result in considerably less accurate reconstruction results. 

At the practical level, the reconstruction is implemented based on the split 
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Bregman approach (with 20 being the maximum number of Bregman iterations 
used in this study). The resulting algorithm alternates between two estimation 
stages (38): first, a sequence of basis pursuit de-noising problems are solved 
independently on a voxel-by-voxel basis, followed by applying a TV filter to a 
total of K discrete images. Such computations are straightforward to acceler- 
ate using a multicore processing, which is another advantage of the proposed 
reconstruction method. 

We believe that the algorithm presented in this paper can be improved in a 
number of ways. First, the square metric used to assess the model fidelity could 
be replaced by a different metric, which would be more specific to the nature of 
Rician noise. Second, the fact that diffusion signals are positive- valued could be 
explicitly incorporated into the reconstruction process in the form of additional 
constraints. Lastly, the bounded variation model could be substituted by an 
alternative model, which could (possibly) provide a better account for the spatial 
regularity of HARDI signals. Exploring the above options constitutes essential 
part of our ongoing research. 

Finally, as the experimental study reported in this paper was comparative in 
its nature, it was not really important what method to use for approximation of 
ODFs. Specifically, the present results have been obtained using Tuch's approx- 
imation 14 . However, more accurate computation of ODFs is possible based on 
the solid angle formulation as detailed in [89, 90|^ It should also be noted that 



the technique proposed in 89 can be applied to multi-shell HARDI data (i.e. 
the data acquired for a range of 6- values). Until recently, collecting such data 
has been deemed impractical due to extremely long acquisition time required. 
We believe, however, that the proposed method for CS-based reconstruction of 
HARDI data has a potential to help multi-shell HARDI develop into a clinically 
relevant tool of diagnostic imaging. 
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